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ABSTRACT 

An  entirely  new  approach  to  the  large-eddy  simulation  (LES)  of  high-speed  compressible 
turbulent  flows  is  presented.  Subgrid  scale  stress  models  are  proposed  that  are  dimensionless 
functions  of  the  computational  mesh  size  times  a  Reynolds  stress  model.  This  allows  a 
DNS  to  go  continuously  to  an  LES  and  then  a  Reynolds-averaged  Navier-Stokes  (RANS) 
computation  as  the  mesh  becomes  successively  more  coarse  or  the  Reynolds  number  becomes 
much  larger.  Here,  the  level  of  discretization  is  parameterized  by  the  non-dimensional  ratio 
of  the  computational  mesh  size  to  the  Kolmogorov  length  scale.  The  Reynolds  stress  model  is 
based  on  a  state-of-the-art  two-equation  model  whose  enhanced  performance  is  documented 
in  detail  in  a  variety  of  benchmark  flows.  It  contains  many  of  the  most  recent  advances 
in  compressible  turbulence  modeling.  Applications  to  the  high-speed  aerodynamic  flows  of 
technological  importance  are  briefly  discussed. 


1.  INTRODUCTION 


The  high-speed  turbulent  flows  of  aerodynamic  importance  exhibit  such  a  wide  range  of 
excited  length  and  time  scales  that  Direct  Numerical  Simulations  (DNS)  are  all  but  impos¬ 
sible  for  the  foreseeable  future.  These  flows  must  invariably  be  computed  based  on  some 
form  of  modeling.  It  has  often  been  said  that  large-eddy  simulations  -  where  the  large  scales 
are  computed  directly  while  the  small  scales  are  modeled  -  are  the  most  promising  such  ap¬ 
proach  (see  Rogallo  and  Moin  1984  for  an  interesting  review).  However,  as  will  be  discussed 
in  this  paper,  there  are  major  problems  with  the  way  that  traditional  large-eddy  simulations 
have  been  conducted.  Furthermore,  there  is  little  experience  with  the  large-eddy  simulation 
of  high-speed  compressible  turbulent  flows  (one  of  the  few  exceptions  is  the  recent  paper 
by  Erlebacher,  Hussaini,  Speziale  and  Zang  1992).  In  addition,  in  practical  wall-bounded 
turbulent  flows,  the  Reynolds  number  can  be  so  high  (in  many  Naval  applications  it  can  be 
as  large  as  109)  that  the  only  feasible  option  is  a  RANS  computation  at  this  time.  Hence, 
there  is  the  need  for  a  combined  LES  and  RANS  capability  for  compressible  turbulent  flows 
where  the  latter  is  of  the  time-dependent  variety.  This  forms  the  motivation  for  the  present 
paper. 

Traditional  large-eddy  simulations  have  typically  been  based  on  the  Smagorinsky  (1963) 
model.  The  Smagorinsky  model  -  along  with  its  extensions  such  as  the  linear  combination 
model  (see  Bardina,  Ferziger  and  Reynolds  1983)  -  have  the  following  two  problems: 

(1)  The  inability  to  respond  to  changes  in  the  local  state  of  the  flow,  causing  the  need  to 
make  ad  hoc  adjustments  in  the  constants. 

(2)  The  generally  poor  correlation  with  DNS  at  lower  turbulence  Reynolds  numbers,  even 
for  simple  benchmark  cases. 

As  far  as  point  (1)  is  concerned,  turbulent  channel  flow,  isotropic  turbulence  and  more 
general,  homogeneously  strained  turbulent  flows  all  require  different  values  of  the  Smagorin¬ 
sky  constant  that  can  differ  by  more  than  a  factor  of  two,  even  when  given  in  its  traditional 
form  of  VCS  where  Cs  is  the  Smagorinsky  constant.  Furthermore,  Van  Driest  wall  damping 
has  been  needed  which  is  empirical  in  nature  and  does  not  apply  to  general  wall-bounded  tur¬ 
bulent  flows  -  particularly  to  those  in  complex  geometries  or  with  flow  separation.  Then,  as 


2 


far  as  point  (2)  is  concerned,  even  for  the  simple  case  of  isotropic  turbulence,  the  Smagorin- 
sky  model  only  correlates  at  the  50%  level  -  an  extremely  poor  result.  It  should  be  noted, 
for  example,  that  the  correlation  between  the  functions  y  —  x  and  y  =  e~x  on  the  interval 
[0,  1]  is  more  than  60%  despite  the  fact  that  they  are  qualitatively  different  functions  (one 
is  monotonically  increasing  while  the  other  is  a  monotonically  decreasing  function)! 

The  only  reason  to  believe  that  the  Smagorinsky  model  is  successful  in  these  cases  is 
probably  due  to  the  fact  that  it  dissipates  enough  energy  to  roughly  account  for  the  cascade 
of  energy  to  the  scales  that  are  left  unresolved.  With  the  desire  to  eliminate  these  problems, 
the  dynamic  subgrid  scale  model  was  recently  developed  (see  Germano,  Piomelli,  Moin  and 
Cabot  1991).  In  the  dynamic  subgrid  scale  model,  a  test  filter  is  introduced  in  addition  to 
the  grid  filter  that  is  used  in  traditional  LES.  A  variable  Smagorinsky  coefficient  is  then 
derived  that  depends  on  the  local  filtered  rate-of-strain  tensor  as  well  as  on  the  resolved 
turbulent  stresses.  The  Smagorinsky  coefficient  then  has  the  capability,  in  principle,  of 
adapting  automatically  to  the  local  state  of  the  flow.  While  the  dynamic  subgrid  scale 
model  does  contain  many  interesting  new  ideas,  it,  unfortunately,  can  further  exacerbate  the 
problem  of  contamination  of  the  large  scales  by  filtering  and  is  not  suitable  for  turbulent 
flows  in  complex  geometries  where  the  effect  of  the  filter  is  never  known  with  certainty  and 
defiltering  is  not  possible  with  any  reliability.  Furthermore,  the  dynamic  subgrid  scale  model 
suffers  from  the  same  deficiency  as  the  older  subgrid  scale  models  since,  in  the  coarse  mesh 
and  infinite- Reynolds-number  limit,  a  state-of-the-art  Reynolds  stress  model  is  not  recovered 
(the  Smagorinsky  model  goes  to  a  badly  calibrated  mixing  length  model  in  the  coarse  mesh 
limit).  In  the  opinion  of  the  author,  future  subgrid  scale  models  must  be  theoretically 
based  on  a  filter  which  does  not  significantly  contaminate  the  large  scales  in  so  far  as  the 
model  calibration  is  concerned  -  with  the  understanding  that  for  complex  turbulent  flows 
one  will  never  know  precisely  the  effect  of  the  filter  and  the  filtered  velocity  must  be  used  to 
approximate  the  large-scale  part  of  the  instaneous  velocity  for  the  calculation  of  turbulence 
statistics.  In  this  manner,  the  issue  of  defiltering  is  completely  avoided  since  it  can  never  be 
done  reliably  anyhow  (defiltering  is  equivalent  to  solving  a  Fredholm  integral  equation  of  the 
first  kind  which  is  mathematically  ill-posed;  see  Hussaini,  Speziale  and  Zang  1990).  Then, 
even  more  importantly,  a  state-of-the-art  Reynolds  stress  model  should  be  recovered  in  the 
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coarse  mesh /infinite- Reyn olds-number  limit.  The  dynamic  subgrid  scale  model  fails  on  both 
counts.  Thus,  new  approaches  must  be  tried. 

An  entirely  new  approach  to  compressible  large-eddy  simulations  will  be  proposed  in 
order  to  bridge  the  gap  between  LES  and  RANS.  This  new  methodology  will  be  based  on 
state-of-the-art  Reynolds  stress  models  that  incorporate  many  of  the  most  recent  develop¬ 
ments  in  compressible  turbulence  modeling  (see  Speziale  1996).  Subgrid  scale  models  will  be 
proposed  that  continuously  go  to  Reynolds  stress  models  in  the  coarse  mesh/infinite  Reynolds 
number  limit.  Furthermore,  they  automatically  vanish  in  the  fine  mesh  limit  where  a  DNS 
is  recovered.  These  points  will  be  discussed  in  detail  in  the  sections  to  follow  and  a  variety 
of  illustrative  calculations  will  be  provided. 
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2.  THE  COMPRESSIBLE  REYNOLDS  STRESS  MODEL 


In  compressible  turbulence,  the  velocity  ut,  pressure  p  and  absolute  temperature  T  are 
solutions  of  the  (cf.  Cebeci  and  Smith  1974) 

Navier-Stokes  Equations 

Q  2 

{pui)  H"  {pUiUj)lj  ~  U3,i)\  V  (^) 

Continuity  Equation 

j-t  +  (puih  =  o  (2) 

Energy  Equation 

—(pCvT)  +  (pCvUiT),i  =  -puiti  +  $  +  (kT,{  )„•  (3) 

where  the  commonly  used  ideal  gas  equation  of  state  is  implemented  in  which 

p  =  pRT. 

Here,  the  viscous  dissipation  function 

2 

$  =  (TijUij  =  +  Uji)u13 

and  p  is  the  mass  density,  R  is  the  ideal  gas  constant,  p  is  the  dynamic  viscosity,  Cv  is  the 
specific  heat  at  constant  volume,  <r,j  =  —\puk,kt>ij  +  +wj,»)  is  the  viscous  stress  tensor 

and  k  is  the  thermal  conductivity.  In  (l)-(3),  the  Einstein  summation  convention  applies  to 
repeated  indices  and  (  •  ),,•  denotes  a  gradient  with  respect  to  the  spatial  coordinates  a;,-. 

Two  alternative  decompositions  for  the  velocity  field  a,,  pressure  p  and  temperature  T 
into  mean  and  fluctuating  parts  can  be  defined: 

m^Ui  +  u'i,  p  =  p  +  p,  T  =  T  +  T'  (4) 

where  an  overbar  represents  a  standard  ensemble  mean,  and 

ui  =  ui  +  u”,  p  =  p  +  p",  T  =  T  +  T"  (5) 

where  a  tilde  represents  a  Favre  or  mass  weighted  average  defined  as 


for  any  flow  variable  T.  A  direct  averaging  of  (l)-(3),  yields  the  Reynolds- averaged  Navier- 
Stokes,  continuity  and  energy  equations  given  by  (cf.  Cebeci  and  Smith  1974  and  Speziale 
1996): 

Reynolds-Averaged  Navier-Stokes  Equations 

q  2  _ 

+  {pUiUj),j  =  -pti  -  +[KUi,j  +  -  ( PTij),j  (6) 

Reynolds- Averaged  Continuity  Equation 

§  +  <J»ih  =  0  (?) 

Reynolds- Averaged  Energy  Equation 

Q 

—(pCvT)  +  (puiCvf)ti  =  -pui,i  -  pu'li 

-p'u'i  i  +  $  +  (kT,),,-  -  Qiti  (8) 

where 

Tij  =  u"u'j,  p  =  pRT  Qi  —  pCvu"T" 

are  the  Favre- averaged  Reynolds  stress  tensor,  the  ensemble  averaged  pressure  and  the  Favre- 
averaged  Reynolds  heat  flux  where  turbulent  fluctuations  in  Cv  have  been  neglected. 

In  high-Reynolds-number  turbulent  flows,  the  molecular  diffusion  terms  are  dominated 
by  the  turbulent  transport  terms  except  in  a  thin  sublayer  near  walls.  If  we  assume  in  this 
region  that  fluctuations  in  the  viscosity,  thermal  conductivity  and  density  can  be  neglected, 
we  can  then  make  the  standard  approximations  (see  Speziale  1996): 


2 

(Tij  =  --flU^kdij  +  p{Uij  +  ujti) 

(9) 

2 

~  -  -puk,k  Si  j  +  p{Ui,]  +  u  jti ) 

Ri  — kT 

(10) 

The  mean  viscous  dissipation  function  is  given  by: 

¥  -  WijUij  +  a'iju'ij 

(11) 

=  +  pe 
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where  e  =  cr'-u'^j/p  is  the  turbulent  dissipation  rate  ($  — *  pe  as  Re  — >  oo).  The  pressure- 
dilatation  correlation,  mass  flux  and  Reynolds  heat  flux  are,  respectively,  modeled  by  (see 
Speziale  1996): 

p'u'ii  =  a2pTlJUijMt  +  azpeMf  (12) 

<  =  (13) 

P  P^P  ^ 

_  p 

(14) 

rrj  e 

where 

C„  =  0.09,  <rp  =  0.5 
Ptj  =  0.9,  o2  =  0.15,  a3  =  0.2. 

and  Mt  =  ( 2K/~fRT )*/2  is  the  turbulent  Mach  number  (see  Speziale  1996,  Speziale  and 
Sarkar  1991  and  Sarkar  1992).  Here,  Ptt  is  the  turbulent  Prandtl  number  and  K  = 
is  the  turbulent  kinetic  energy;  7  =  CPICV  is  the  ratio  of  mean  specific  heats  at  constant 
pressure  and  constant  volume. 

The  Reynolds  stress  is  represented  by  the  explicit  algebraic  stress  model  (see  Speziale 
1996): 

Tij  =  \KSij  -  -  ^SkkStj)  -  a*2^-(Sikukj 

J  1  '  (15) 

~\~Sjkivki')  T  (sikSkj  ” Ski Skl 


where 


for  i  =  1, 2, 3  and 


a-  —  a i 


3  —  2t/2  +  6£2 


a.  =  (|  -  C2)  g,  aj  =  1  (j  -  C2)  (2  -  Ct)g 2 

a3  =  (g  —  £2)  (2  —  CS)#2,  g  =  (2^3  +  —  -  l) 

V  =  i^(S,A)1/2,  f  =  “(* (U) 
z  ol\  £  ol\  e 

This  explicit  algebraic  stress  model  -  which  is  in  the  form  of  an  anisotropic  eddy  viscosity 
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model  with  strain-dependent  coefficients  -  is  obtained  from  the  explicit  solution  of  the  equi¬ 
librium  form  of  the  modeled  Reynolds  stress  transport  equation  for  incompressible  turbulent 
flows  (see  Gatski  and  Speziale  1993).  It  is  then  extended  to  compressible  flows.  In  (17),  V /e 
is  set  to  one  and  C\  —  C4  are  constants  that  assume  the  values  of  5.20,  0.36,  1.25  and  0.40, 
respectively.  Here,  again  the  turbulent  kinetic  energy  K  =  |r,;.  For  general  turbulent  flows, 
a*  is  represented  by  the  regularized  expressions 


(1  +  2^(1  +  6^5)  +  h2 

ai  (l  +  2£2)(l  +  2£2  +r12  +  6plrie)  1 

.  (l  +  2C2)(l  +  V)  +  h2 

(1  +  2£2)(1  +  2£2  +  /3,-7/6)a* 

for  i  =  2, 3  where 

Px  =  7.0,  /32  =  6.3,  ft  =  4.0. 


(18) 

(19) 


This  is  obtained  by  a  Pade’  approximation  which  builds  in  some  limited  agreement  with 
Rapid  Distortion  Theory  (RDT)  (see  Speziale  and  Xu  1996). 

The  turbulent  kinetic  energy  is  a  solution  of  the  transport  equation: 


Q 

—(pK)  +  (puiK)ti  =  —pTijUij  —  pe 


+P'K,i  ~  UiP,i  + 

+  f(V+-Kl 

LV  Ck J  J,,- 


(20) 


where 


Pt  —  C»P 


K 2 


is  the  eddy  viscosity  and  C ^  «  0.09  is  taken  to  be  a  constant;  <Tk  is  a  constant  equal  to  one. 
The  dissipation  rate  equation  takes  the  compressible  modeled  form  (see  Speziale  1996): 


d_ 

dt 


(pe)  +  (pui£),i  =  1  p ~~ T{j  (iiij  -  ^uk,k$ij^ 


■  —p£Ui  i  +  CtfRt1/2—  —  Ce2p-jg 


(21) 


+ 
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where  ae  is  a  constant  that  assumes  the  value  of  approximately  1.3  and  Ce\  =  1.44,  Ce 2  =  1.83 
and  Ce3  =  0.001;  Rt  =  K2 I've  is  the  turbulence  Reynolds  number.  This  dissipation  rate 
model  can  be  directly  integrated  to  a  solid  boundary  with  no  wall  damping  (see  Speziale 
and  Abid  1995).  Only  the  singularity  need  be  removed  in  the  destruction  term  where  a 
simple  modification  in  CE2  is  made.  The  vortex  stretching  term  in  the  expression  containing 
C£ 3  leads  to  a  better  behaved  model  without  compromising  the  predictions  in  benchmark 
equilibrium  turbulent  flows.  For  example,  it  removes  the  singularity  in  plane  stagnation  point 
turbulent  flows  (see  Abid  and  Speziale  1996).  Furthermore,  it  allows  a  universal  equilibrium 
value  of  V/e  —  1  to  be  used  in  the  explicit  algebraic  stress  models  with  good  agreement  for 
both  homogeneous  shear  flow  and  the  log-layer.  Setting  V/e  to  a  constant  makes  the  model 
more  robust.  The  Sarkar  compressible  dissipation  model  (see  Sarkar  et  al.  1991  and  Zeman 
1990)  was  not  implemented  because  it  is  either  negligible  or  does  not  do  well  in  the  log-layer 
and  it  is  our  purpose  to  develop  a  model  suitable  for  supersonic  wall-bounded  turbulent 
flows. 
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3.  A  NEW  APPROACH  TO  LARGE-EDDY  SIMULATIONS 


In  the  large-eddy  simulation  of  compressible  turbulent  flows  the  filtered  equations  of 
motion  are  given  by: 

Filtered  Navier-Stokes  Equations 

3  2 

—  (pzii)  +  (pUiUj)j  —  ~Pj  —  ~~  {prij),j  (22) 


Filtered  Continuity  Equation 

Yt  +  (?“•■),•■  = 0  (23) 

Filtered  Energy  Equation 

3 

T-(pc.r)  +  (rm.c.f),  =  -pii,,  -  ?«,)s 

+  5s  +  (JSTi),.-  -  Ql  (24) 

_  p 

where  r-,  (w")5,  (p'u(- ,)s,  $  and  are,  respectively,  the  subgrid  scale  Reynolds  stress  ten¬ 
sor,  the  subgrid  scale  mass  flux,  the  subgrid  scale  pressure-dilatation  correlation,  the  subgrid 
scale  dissipation  function  and  the  subgrid  sea  le  Reynolds  heat  flux.  The  approximations 
discussed  in  the  previous  section  for  the  diffusion  terms  in  (9)— (10)  were  made.  It  should  be 
noted  that  the  filtered  equations  of  motion  are  identical  in  form  to  the  Reynolds-averaged 
equations  of  motion  discussed  earlier.  This  will  allow  us  to  develop  a  combined  LES  and 
time- dependent  RANS  capability.  In  (22)-(24)  an  overbar  represents  a  standard  filter  and 
a  tilde  represents  a  Favre  filter.  For  any  flow  variable  J-  they  are  defined  as  follows: 


T=  [  G(x  —  x*,  A)Jr(x*)d3x* 

J  D 


(25) 


where  G  is  a  filter  function,  A  is  the  computational  mesh  size  and 


A  compact  filter  of  bounded  support,  such  as  the  box  filter,  can  be  implemented  in  complex 
geometries  where 


G(x  —  x*,  A) 


1/8A3,  \x(  —  a:,*|  <  A 
0,  |x,-  —  ®,*|  >  A 


(26) 
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for  i  =  1,2,3  (cf.  Leonard  1974  and  Ferziger  1976). 

Models  for  the  subgrid  scale  stress  tensor  rfj  are  proposed  that  are  of  the  form: 

r?  =  [1  -  exP(-^A/iK)]"ry  (27) 

where  r,j  is  a  Reynolds  stress  model,  A  is  the  computational  mesh  size,  Lk  is  the  Kolmogorov 
length  scale,  and  /3  and  n  are  constants.  In  the  Reynolds  stress  model,  the  ensemble  averages 
are  taken  to  be  filtered  quantities  (in  complex  geometries  we  never  know  the  effect  of  the 
filter  anyhow).  In  the  limit  as  A /Lk  — >  0,  all  relevant  scales  are  resolved  and  we  have 
a  direct  simulation  where  r,j  =  0;  as  A/Lk  —■ ►  oo  and  the  mesh  becomes  coarse  or  the 
Reynolds  number  becomes  extremely  large,  we  recover  a  Reynolds  stress  model  and  a  RANS 
computation.  In  between  these  two  limits,  we  have  an  LES  or  a  VLES  (the  latter  denotes  a 
large-eddy  simulation  where  the  preponderance  of  the  turbulent  kinetic  energy  is  unresolved). 
Here,  it  should  be  noted  that  in  the  simulation  of  turbulence,  the  computational  mesh  is 
fine  (or  coarse)  depending  on  whether  A  is  small  (or  large)  compared  to  the  Kolmogorov 
length  scale  Lk  =  i/^/e1/4  where  V  =  ~p/p.  This  automatically  brings  in  a  dependence  on 
the  turbulence  Reynolds  number  Rt  ( Rt  =  K2  / Ve)  since  Lk  =  R^3^4K3^2 /e.  An  estimate  of 
the  Kolmogorov  length  scale  is  provided  by  the  modeled  transport  equation  for  e  discussed 
earlier.  In  order  to  estimate  the  the  Kolmogorov  length  scale  to  within  10%  it  is  only 
required  that  the  dissipation  rate  be  estimated  to  within  50%  since  the  dissipation  rate  is 
raised  to  the  1/4  power.  Hence,  it  is  quite  reasonable  to  expect  that  an  acceptable  estimate 
of  the  Kolmogorov  length  scale  can  be  obtained  (this  is  a  crucial  element  of  this  proposed 
approach).  The  Reynolds  stress  model  in  (27)  is  represented  by  the  expressions  discussed  in 
Section  2.  It  was  recently  shown  by  Speziale  and  Abid  1995  that  this  type  of  two-equation 
model  can  be  integrated  directly  to  a  solid  boundary  with  no  wall  damping  functions;  only 
the  singularity  in  the  e  -  transport  equation  needs  to  be  removed. 

Initially  we  have  chosen 

ft  «  0.001,  n  =  1. 

Hence,  for  A/Lk  <  10  -  which  is  a  value  encountered  in  many  practical  LES  -  the  grid 
function  in  (27)  can  be  approximated  by  a  power  law.  Interestingly  enough,  Woodruff  and 
Hussaini  ( Private  Communication )  have  arrived  at  a  power  law  representation  for  the  grid 
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function  in  (27)  via  Renormalization  Group  techniques.  This  shows  that  there  is  some 
theoretical  basis  for  this  idea. 

The  same  modeling  technique  is  implemented  for  the  subgrid  scale  mass  flux,  subgrid 
scale  pres  sure- dilatation  correlation,  subgrid  scale  dissipation  function  and  subgrid  scale 
Reynolds  heat  flux.  More  specifically,  subgrid  scale  models  are  implemented  which  are  of 
the  form 


W)S  =  [1  -  e*?(-/?A/IK)]X 

(28) 

(P'<.r  =  f1  -  ea:p(-/3A/r^)]V<i 

(29) 

=  [1  -  exp{-/3A/LK)]n$ 

(30) 

Qf  =  [  1  -  exp(-(3A /LK)]nQi 

(31) 

where  the  mass  flux,  pressure-dilatation  correlation,  dissipation  function  and  Reynolds  heat 
flux  on  the  right-hand-side  of  (28)-(31)  are  modeled  in  the  Reynolds-averaged  forms  pro¬ 
vided  in  Section  2.  However,  an  overbar  or  tilde  is  taken  to  be  a  standard  or  Favre  filter, 
respectively.  For  the  practical  LES  of  high  Reynolds  number  turbulent  flows,  the  subgrid 
scale  mass  flux  and  subgrid  scale  pressure-dilatation  correlation  can  be  neglected  (see  Er- 
lebacher,  Hussaini,  Speziale  and  Zang  1992).  Furthermore,  the  mean  dissipation  function 
can  be  approximated  by  $  «  pe  which  leads  to  a  considerable  simplification. 

Some  comments  are  needed  concerning  the  choice  of  a  filter  in  this  new  approach  to 
large-eddy  simulations.  We  want  a  filter  that  yields  the  minimum  contamination  of  the  large 
scales.  The  reason  for  this  is  clear;  defiltering  must  be  avoided  since  it  constitutes  an  ill-posed 
mathematical  problem  as  stated  earlier.  The  main  purpose  of  practical  LES  is  to  predict 
the  Reynolds- averaged  fields.  In  order  to  do  so,  the  filtered  velocity  must  invariably  be  used 
to  estimate  the  instaneous  velocity  which  then  yields  the  Reynolds- averaged  fields  through 
appropriate  ensemble  averages.  Here,  the  large  scales  make  the  dominant  contribution  to 
the  most  pertinent  fields  such  as  the  turbulent  kinetic  energy.  A  minimum  contamination 
of  the  large  scales  can  be  accomplished  with,  of  the  order  of,  a  1283  computational  mesh 
using  a  filter  with  a  compact  support  -  such  as  the  box  filter  -  which  has  a  small  filter  width 
of,  say,  two  mesh  points.  Some  of  the  previously  conducted  coarse  grid  LES  (which  has 
typically  had  no  more  than  323  mesh  points)  must  be  avoided  wherein  the  filter  width  has, 
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at  times,  been  as  much  as  25%  of  the  computational  domain,  significantly  contaminating  the 
large  scales.  Besides,  recent  increases  in  computational  capacity  have  begun  to  make  1283 
computations  much  more  feasible  for  engineering  calculations  (a  small  compromise  to  1003 
computations  can  always  be  made).  In  addition,  it  should  be  noted  that  practical  LES  -  in 
complex  geometries  -  will  require  the  use  of  finite  difference  techniques  with  a  compact  filter 
(these  should  be  based  on  fourth-order  accurate  finite  difference  schemes).  Some  illustrative 
calculations  will  be  provided  in  the  next  Section  to  demonstrate  the  potential  of  this  new 
approach  for  practical  engineering  calculations. 
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4.  ILLUSTRATIVE  CALCULATIONS 


First,  we  will  present  some  results  for  incompressible  flows  to  demonstrate  that  the  model 
behaves  correctly  in  the  limit  of  small  Mach  numbers.  As  mentioned  earlier,  this  two-equation 
model  can  be  integrated  directly  to  a  solid  boundary  with  no  wall  damping.  Results  were 
obtained  for  the  incompressible  flat  plate  boundary  layer  at  zero  pressure  gradient.  In 
Figure  1,  the  skin  friction  coefficient  obtained  from  this  model  -  plotted  as  function  of  the 
Reynolds  number  based  on  the  momentum  thickness,  Re  -  is  compared  with  experimental 
data  (see  Wieghardt  and  Tillman  1951)  and  with  results  obtained  from  the  K  —  e  model  with 
wall  damping  (see  Speziale,  Abid  and  Anderson  1992).  Clearly,  the  results  are  extremely 
good.  No  wall  damping  functions  are  needed;  the  singularity  which  occurs  at  the  wall  only 
needs  to  be  removed  from  the  dissipation  rate  equation  (see  Speziale  and  Abid  1995).  The 
fact  that  excellent  predictions  for  the  Reynolds  stresses  are  obtained  starting  in  the  log- 
layer  is  demonstrated  in  Figure  2  where  the  turbulence  intensities  are  depicted.  A  Reynolds 
stress  model  cannot  be  expected  to  accurately  describe  the  turbulence  structure  immediately 
adjacent  to  a  wall  since  it  exhibits  such  a  wide  range  of  scales  and  diverse  turbulence  physics. 
We  only  need  to  have  Reynolds  stress  models  that  can  operationally  be  integrated  to  a 
wall  without  compromising  their  predictive  capabilities  away  from  the  wall.  This  can  be 
accomplished  with  the  model  discussed  herein. 

In  Figure  3,  the  prediction  of  this  two-equation  model  for  the  mean  velocity  profile  in 
rotating  channel  flow  is  compared  with  the  experimental  data  of  Johnston,  Halleen  and 
Lezius  (1972)  for  a  rotation  number  Ro  =  0.068.  It  is  clear  from  these  results  that  the  model 
correctly  predicts  that  the  mean  velocity  profile  is  asymmetric  in  line  with  the  experimental 
results  -  an  effect  that  arises  from  Coriolis  forces.  In  contrast  to  these  results,  the  standard 
K  —  e  model  incorrectly  predicts  a  symmetric  mean  velocity  profile  identical  to  that  obtained 
in  an  inertial  frame  (the  standard  K—e  model  is  oblivious  to  rotations  of  the  reference  frame). 
As  demonstrated  by  Gatski  and  Speziale  (1993),  the  results  obtained  in  Figure  3  with  this 
two-equation  model  are  virtually  as  good  as  those  obtained  from  a  full  second-order  closure. 
This  is  due  to  the  fact  that  a  representation  is  used  for  the  Reynolds  stress  tensor  that  is  a 
two-equation  model  which  is  formally  derived  from  a  second-order  closure  in  the  equilibrium 
limit.  It  is  now  clear  that  previous  claims  that  two-equation  models  cannot  systematically 


14 


account  for  rotational  effects  were  erroneous. 

A  few  examples  will  now  be  presented  that  illustrate  the  enhanced  predictions  that  are 
obtained  for  turbulent  flows  exhibiting  effects  arising  from  normal  Reynolds  stress  differences. 
Here,  we  will  show  results  obtained  from  the  nonlinear  K  —  e  model  of  Speziale  (1987)  (also 
see  Yoshizawa  1984).  For  turbulent  shear  flows  that  are  predominantly  unidirectional,  with 
secondary  flows  or  recirculation  zones  driven  by  small  normal  Reynolds  stress  differences,  a 
quadratic  approximation  to  the  anisotropic  eddy  viscosity  model  described  herein  collapses 
to  the  nonlinear  K  —  e  model  (see  Gatski  and  Speziale  1993).  In  Figure  4,  it  is  demonstrated 
that  the  nonlinear  K  —  e  model  predicts  an  eight-vortex  secondary  flow,  in  a  square  duct,  in 
line  with  experimental  observations;  on  the  other  hand,  the  standard  K—e  model  erroneously 
predicts  that  there  is  no  secondary  flow.  In  order  to  be  able  to  predict  secondary  flows  in  non- 
circular  ducts,  the  axial  mean  velocity  uz  must  give  rise  to  a  non-zero  normal  Reynolds  stress 
difference  ryy  —  txx.  This  requires  an  anisotropic  eddy  viscosity  (any  isotropic  eddy  viscosity, 
including  that  used  in  the  standard  K  —  e  model,  yields  a  vanishing  normal  Reynolds  stress 
difference  which  makes  it  impossible  to  describe  these  secondary  flows). 

In  Figure  5,  results  obtained  from  the  nonlinear  K  —  e  model  are  compared  with  the 
experimental  data  of  Kim,  Kline  and  Johnston  (1980)  and  Eaton  and  Johnston  (1980)  for 
turbulent  flow  past  a  backward  facing  step.  It  is  clear  that  these  results  are  excellent: 
reattachment  is  predicted  at  x/H  ~  7.0  in  close  agreement  with  the  experimental  data 
(see  Figure  5(a)).  The  Reynolds  shear  stresses  are  also  well  represented  by  the  model  (see 
Figure  5(b)).  In  contrast  to  these  results,  the  standard  K  —  e  model  predicts  reattachment  at 
x/H  ps  6.25  -  an  11%  underprediction.  This  error  predominantly  results  from  the  inaccurate 
prediction  of  normal  Reynolds  stress  anisotropies  in  the  recirculation  zone  as  discussed  by 
Speziale  and  Ngo  (1988). 

The  basic  two-equation  model  discussed  herein  has  recently  been  shown  to  be  feasible  for 
complex  aerodynamic  computations  by  Abid,  Morrison,  Gatski  and  Speziale  (1996)  in  the 
compressible  regime  (the  transonic  flow  regime  to  be  specific).  The  ONERA  M6  wing  at  a 
Mach  number  of  M0 0  =  0.8447  with  an  angle  of  attack  of  a  =  5.06°  and  a  Reynolds  number  of 
Re  =  11.7  x  106  has  been  considered.  At  this  Mach  number  the  turbulence  statistics  behave 
quasi-incompressibly  (compressibility  effects  are  primarily  felt  through  changes  in  the  mean 
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density  since  the  Morkovin  hypothesis  applies).  In  Figure  6,  a  comparison  is  made  of  the 
computed  surface  pressure  distributions  with  experiments  (see  Schmitt  and  Charpin  1979)  at 
four  different  spanwise  locations.  As  can  be  seen,  the  results  are  quite  good.  They  are  shown 
primarily  to  establish  the  feasibility  of  this  new  two-equation  model  for  the  calculation  of 
complex  compressible  flows. 

Some  preliminary  LES  results  will  now  be  provided  for  the  developing  incompressible 
turbulent  boundary  layer  at  zero  pressure  gradient  —  integrated  through  transition.  These 
computations  were  conducted  by  H.  Fasel  and  his  group  at  the  University  of  Arizona  using 
an  empirically  based  ramp  function  —  that  depends  explicitly  on  the  momentum  thickness 
Reynolds  number  and  the  mesh  size  with  a  simple  eddy  viscosity  model  —  as  a  preliminary 
test  of  the  ideas  embodied  in  (27).  In  Figure  7,  the  spanwise  vorticity  obtained  from  the  LES 
is  shown  which  compares  favorably  with  the  corresponding  results  obtained  from  DNS.  It  is 
clear  that  the  subgrid  scale  model  allows  the  LES  to  pick  up  the  pertinent  flow  structures 
and  to  be  integrated  through  transition  (laminar  -  turbulent  flow).  The  ramp  function, 
which  forms  a  central  part  of  (27),  allows  the  eddy  viscosity  to  gradually  turn  on  as  the  flow 
becomes  turbulent.  In  this  regard,  the  corresponding  eddy  viscosity  is  displayed  in  Figure 
8.  These  calculations  -  which  have  yielded  encouraging  results  -  are  being  repeated  for  the 
compressible  turbulent  boundary  layer. 

Unlike  in  many  existing  dissipation  rate  models,  the  dilatational  part  of  the  production 
terms  in  (21)  are  consistent  with  rapid  distortion  theory  for  compressible  isotropic  turbulence. 
For  the  problem  of  a  rapid  compression  or  expansion  of  isotropic  turbulence,  the  mean 
velocity  gradient  tensor  is  given  by 
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where  T  is  the  expansion/compression  rate  which  is  constant.  The  two-equation  model 
provided  herein  reduces,  in  this  case,  approximately  to  the  simple  system  of  coupled  ODE’s: 


K  =  -\vK  (33) 

e  =  -^re  (34) 
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for  \T\Ko/eo  1.  The  short-time  solution  to  Eqs.  (33)  -  (34)  is  given  by: 


K  =  Kq  exp 

(35) 

e  =  £0  exp 

(36) 

A  =  A0  exp 

(37) 

where  A  =  K:i/2 /e  is  the  integral  length  scale.  This  solution  is  essentially  identical  to  the 
results  obtained  by  Reynolds  (1987)  based  on  Rapid  Distortion  Theory  (RDT).  In  contrast 
to  these  results,  a  variable  density  extension  of  the  commonly  used  second-order  closures 
erroneously  predicts  that 

A  =  Aoexp(-0.04rt)  (38) 

(see  Reynolds  1987).  According  to  (38),  the  integral  length  scale  will  decrease  under  an 
expansion  (T  >  0)  and  increase  under  a  compression  (1?  <  0)  -  results  that  are  clearly  in 
error  as  first  pointed  out  by  Reynolds  (1987).  It  is  believed  that  this  correction  will  lead  to 
a  better  behaved  model  in  regions  of  strong  straining  which  can  occur  when  shock  waves  are 
present. 

Interesting  enough,  good  results  are  obtained  for  the  compressible  flat  plate  boundary  layer 
with  this  two-equation  model  for  Mach  numbers  Moo  as  large  as  8  and  wall  temperature  ratios 
as  low  as  0.3  with  no  explicit  dilatation  terms,  as  first  calculated  by  Zhang  et  al.  (1993) 
(see  Figure  9).  In  fact,  it  must  be  noted  that  the  compressible  dissipation  model  of  Sarkar 
causes  a  degradation  of  the  skin  friction  predictions;  the  model  is  not  formally  correct  in 
the  log-layer  of  high-speed  compressible  boundary  layers.  This  was  probably  first  pointed 
out  by  P.  G.  Huang  (see  Huang  et  al.  1994).  That  is  why  we  have  chosen  to  exclude  the 
compressible  or  dilatational  dissipation  from  the  model.  It  is  not  needed  to  get  good  results 
in  the  compressible  turbulent  boundary  layer  and  we  are  interested  in  describing  high-speed 
compressible  wall-bounded  flows. 
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5.  CONCLUSION 


A  new  combined  LES  and  time-dependent  RANS  capability  for  the  computation  of  the 
high-speed  compressible  turbulent  flows  of  technological  importance  has  been  proposed.  In 
this  new  methodology,  subgrid  scale  models  go  continuously  to  Reynolds  stress  models  in 
the  coarse  mesh/infinite  Reynolds  number  limit.  Here,  it  is  believed  that  in  complex  wall- 
bounded  turbulent  flows  -  especially  with  flow  separation  where  wall  functions  cannot  be 
used  -  the  best  that  one  can  do,  at  this  time,  for  the  extremely  high  Reynolds  numbers  en¬ 
countered  in  many  technological  applications  (such  as  some  naval  flows  where  Re  ~  0(1O9)) 
is  a  RANS  computation  since  the  crucial  wall  layer  cannot  be  resolved.  At  more  moderate 
Reynolds  numbers,  LES  is  possible  (of  course  at  extremely  low  turbulence  Reynolds  num¬ 
bers  it  is  possible  to  conduct  a  DNS).  The  RANS  model  that  was  presented  in  this  study  is 
a  two-equation  model  that  contains  some  of  the  most  recent  developments  in  compressible 
turbulence  modeling.  In  the  incompressible  limit  it  collapses  to  an  explicit  algebraic  stress 
model  which  is  a  two-equation  model  that  is  consistent  with  a  state-of-the-art  second-order 
closure  model  in  the  limit  of  homogeneous  turbulence  in  equilibrium.  A  mild  non-equilibrium 
extension  is  then  built  in  via  a  Pade'  approximation  that  has  some  limited  agreement  with 
Rapid  Distortion  Theory  (RDT).  The  performance  of  this  two-equation  model  was  docu¬ 
mented  in  a  variety  of  benchmark  turbulent  flows.  This  performance  was  surprisingly  good. 

A  few  words  are  warranted  concerning  supersonic  turbulent  flows  -  particularly  in  wall- 
bounded  geometries.  Since  there  is  considerable  heating  in  such  flows,  the  effective  Reynolds 
number  (through  a  rise  in  the  kinematic  viscosity)  is  not  extremely  large  even  in  many  of 
the  high-speed  compressible  flows  of  technological  importance.  Thus,  the  prospects  for  being 
able  to  conduct  a  large-eddy  simulation  in  many  practical  flow  situations  are  fairly  good  for 
compressible  turbulence.  This  would  include  the  flow  around  high-speed  aircraft  and  missiles 
to  name  just  a  few  applications.  For  those  cases  where  LES  is  not  feasible,  a  time-dependent 
RANS  can  be  always  be  conducted  which  is  expected  to  be  far  superior  to  traditional  steady 
RANS.  The  methodology  discussed  in  this  paper  has  the  potential  to  make  a  significant 
impact  on  such  problems  in  each  of  these  dual  areas.  Preliminary  results  are  quite  promising 
and  further  tests  are  currently  underway. 
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Figure  1.  Comparison  of  the  predictions  of  the  proposed  two-equation  model  for  skin  friction 
with  experimental  data  (Wieghardt  and  Tillman  1951)  for  the  incompressible  flat  plate 
turbulent  boundary  layer. 


Figure  2.  Incompressible  turbulent  flat  plate  boundary  layer:  Comparison  of  the  predictions 

of  the  ( - )  proposed  two-equation  model  and  the  (-  -  -)  Speziale,  Abid  and  Anderson 

(1992)  model  with  (  A  )  experimental  data  (Wieghardt  and  Tillman  1951).  (a)  Streamwise 
turbulence  intensity,  (b)  turbulence  intensity  normal  to  wall  and  (c)  turbulent  shear  stress. 


□  Experimental  Data 


—  New  Explicit  ASM 


y/H 


Figure  3.  Comparison  of  the  mean  velocity  profile  in  incompressible  rotating  channel  flow 
predicted  by  the  proposed  two-equation  model  with  the  experimental  data  of  Johnston, 
Halleen  and  Lezius  (1972). 


Figure  4.  Incompressible  turbulent  secondary  flow  in  a  rectangular  duct:  (a)  experiments, 
(b)  standard  K  —  e  model,  and  (c)  nonlinear  K  —  e  model  of  Speziale  (1987). 


(o  experimental  mean  reattachment  point) 

(a) 


—  Computations  with  anisotropic  eddy  viscosity; 
o  Experiments  of  Kim  etal,  1980;  Eaton  &  Johnston,  1980 


Figure  5.  Incompressible  turbulent  flow  past  a  backward  facing  step;  comparison  of  the 
predictions  of  the  nonlinear  K-e  model  of  Speziale  (1987)  with  experiments,  (a)  Streamlines 
and  (b)  turbulent  shear  stress  profiles. 


Figure  6.  Compressible  surface  pressure  distributions  for  the  ONERA  M6  wing  at  four 
different  spanwise  locations  predicted  by  the  proposed  two-equation  model  (  ).  Experiments 

of  Schmitt  and  Charpin  (1979)  denoted  by  (  o  ). 


Figure  7.  Plot  of  spanwise  vorticity  in  the  developing  incompressible  turbulent  boundary 
layer  obtained  from  LES  (computations  done  by  H.  Fasel  and  co-workers  at  the  University 
of  Arizona). 


120  140  160  180  200 


Figure  8.  Plot  of  the  eddy  viscosity  in  the  developing  incompressible  turbulent  boundary 
layer  obtained  from  LES  (computations  done  by  H.  Fasel  and  co-workers  at  the  University 
of  Arizona). 


(a) 


Figure  9.  Comparison  of  the  predictions  of  the  new  two-equation  model  as  first  computed 
by  Zhang  et  al.  (1993)  for  the  compressible  flat  plate  boundary  layer  with  the  experimental 
data  of  Kussoy  and  Horstman  (1991)  for  =  8.18  and  TJT„  =  0.3:  (a)  Mean  velocity 
and  (b)  mean  temperature. 
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